Trabajo Final - Estadistica Espacial

Instalamos los paquetes necesarios

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(sf)

Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1

library(gstat)
library(geoR)

'RandomFieldsUtils' will NOT use OMP
'RandomFields' will NOT use OMP

--------------------------------------------------------------
 Analysis of Geostatistical Data
 For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
 geoR version 1.8-1 (built on 2020-02-08) is now loaded
--------------------------------------------------------------

library(spdep)

Loading required package: sp
Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with:
`install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source')`

library(DataExplorer)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(psych)
library(ggcorrplot)

Loading required package: ggplot2

Attaching package: ‘ggplot2’

The following objects are masked from ‘package:psych’:

    %+%, alpha

library(biogeo)

Attaching package: ‘biogeo’

The following object is masked from ‘package:spData’:

    world

library(ggplot2)
library(leaflet)

Cargamos los datasets

estaciones <- read.csv("data_smn/preprocessed/estaciones_smn_v2.csv")
horarios <- read.csv("data_smn/preprocessed/datohorario20210420.csv")

Observamos como esta compuesto el dataset

glimpse(estaciones)

Rows: 123
Columns: 9
$ NOMBRE           <chr> "BASE BELGRANO II", "BASE CARLINI (EX JUBANY)", "BASE ESPERANZA", "BASE MARAMBIO", "BASE ORC…
$ PROVINCIA        <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "BUENOS AIRES"…
$ LATITUD_GRADOS   <int> 77, 62, 63, 64, 60, 68, 36, 38, 37, 36, 34, 38, 37, 36, 34, 34, 34, 34, 36, 37, 34, 34, 34, …
$ LATITUD_MINUTOS  <int> 52, 14, 24, 14, 45, 8, 50, 44, 43, 12, 32, 1, 26, 21, 36, 49, 33, 58, 2, 56, 33, 41, 40, 27,…
$ LONGITUD_GRADOS  <int> 34, 58, 56, 56, 44, 67, 59, 62, 59, 61, 58, 61, 61, 57, 58, 58, 60, 57, 59, 57, 58, 58, 58, …
$ LONGITUD_MINUTOS <int> 34, 38, 59, 43, 43, 8, 53, 10, 47, 4, 40, 20, 53, 44, 36, 32, 55, 54, 8, 35, 49, 44, 38, 53,…
$ ALTURA           <int> 256, 11, 24, 198, 12, 7, 147, 83, 207, 94, 26, 247, 233, 9, 12, 20, 81, 23, 36, 21, 32, 36, …
$ NRO              <int> 89034, 89053, 88963, 89055, 88968, 89066, 87641, 87750, 87649, 87640, 87570, 87683, 87637, 8…
$ NroOACI          <chr> "SAYB", "SAYJ", "SAYE", "SAWB", "SAYO", "SAYS", "SAZA", "SAZB", "SAZJ", "SAZI", "SADO", "SAB…

glimpse(horarios)

Rows: 2,041
Columns: 8
$ FECHA  <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20…
$ HORA   <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 0, 1, 2, 3, 4, 5…
$ TEMP   <dbl> 20.7, 19.8, 19.5, 19.2, 19.2, 18.9, 19.0, 18.8, 19.0, 19.1, 19.9, 20.3, 21.8, 22.9, 23.2, 23.7, 23.2, …
$ HUM    <int> 70, 79, 79, 81, 81, 84, 84, 84, 84, 84, 81, 82, 68, 72, 68, 66, 71, 68, 75, 72, 73, 68, 68, 68, 90, 88…
$ PNM    <dbl> 1020.8, 1020.9, 1020.9, 1020.4, 1020.0, 1020.3, 1020.7, 1021.1, 1021.6, 1022.2, 1022.4, 1021.9, 1021.7…
$ DD     <int> 70, 70, 70, 70, 70, 50, 70, 50, 50, 50, 50, 50, 50, 70, 70, 70, 70, 70, 90, 90, 90, 90, 70, 70, 50, 50…
$ FF     <chr> "28", "28", "26", "26", "22", "22", "17", "17", "15", "19", "20", "24", "20", "19", "20", "17", "15", …
$ NOMBRE <chr> "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPARQUE AERO", "AEROPA…

Observamos que la variable que representa la velocidad del viento en km/h es de tipo Char, por lo cual debemos convertila en int.

horarios$FF <- as.numeric(horarios$FF)

NAs introduced by coercion

summary(estaciones$ALTURA)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    3.0    53.5   145.0   331.5   455.0  3459.0 

summary(horarios$HORA)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    7.00   12.00   12.16   18.00   23.00 

summary(horarios$TEMP)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 -20.00   16.10   20.20   19.09   23.80   33.60       1 

summary(horarios$HUM)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  18.00   59.00   72.00   71.18   85.00  100.00       4 

summary(horarios$PNM)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  960.8  1011.0  1015.1  1075.3  1019.2  3133.0      23 

summary(horarios$DD)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    0.0    50.0    70.0   118.9   160.0   990.0       1 

summary(horarios$FF)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   0.00    9.00   15.00   16.06   22.00  107.00       1 

hist(estaciones$ALTURA, main = "Histograma de la altura", xlab = "Altura")

hist(horarios$HORA, main = "Histograma de horarios", xlab = "Horarios")

hist(horarios$TEMP, main = "Histograma de temperatura", xlab = "Temperatura")

hist(horarios$HUM, main = "Histograma de humedad", xlab = "Humedad")

hist(horarios$PNM, main = "Histograma de presion atmosferica", xlab = "Presion Atmosferica")

hist(horarios$DD, main = "Histograma de la direccion del viento", xlab = "Direccion del viento")

hist(horarios$FF, main = "Histograma de la velocidad del viento", xlab = "Velocidad del viento")

Creamos unos boxplot para visualizar la distribucionde datos que potencialmente nos interecen para proceder con el analisis

boxplot(estaciones$ALTURA, main = "Boxplot de la altura", xlab = "Altura")

boxplot(horarios$HORA, main = "Boxplot de horarios", xlab = "Horarios")

boxplot(horarios$TEMP, main = "Boxplot de temperatura", xlab = "Temperatura")

boxplot(horarios$HUM, main = "Boxplot de humedad", xlab = "Humedad")

boxplot(horarios$PNM, main = "Boxplot de presion atmosferica", xlab = "Presion Atmosferica")

boxplot(horarios$DD, main = "Boxplot de la direccion del viento", xlab = "Direccion del viento")

boxplot(horarios$FF, main = "Boxplot de la velocidad del viento", xlab = "Velocidad del viento")

#plot_boxplot(horarios, by = "HORA")
#plot_boxplot(horarios, by = "TEMP")

Ahora, veamos que tan correlacionadas estan estas variables

corr <- cor(horarios[, c(2,3,4,5,6,7)], use = "complete.obs")

ggcorrplot(corr, type = "lower", outline.col = "black",
 lab=TRUE,
 ggtheme = ggplot2::theme_gray,
 colors = c("#6D9EC1", "white", "#E46726"))

Hay datos nulos en estos datasets?

print("Nulos en horarios")

[1] "Nulos en horarios"

which(is.na(horarios))

 [1]  4853  6710  6894  7999  8000  8401  8935  9137  9138  9139  9140  9141  9142  9143  9144  9556  9557  9558  9559
[20]  9560  9561  9562  9563  9564  9565  9566  9567  9605 10976 13017

print("Nulos en estaciones")

[1] "Nulos en estaciones"

which(is.na(estaciones))

integer(0)

Borramos los nulos del dataset

estaciones = na.omit(estaciones)
horarios = na.omit(horarios)

Analizamos la simetria de la variable que representa la velocidad del vientoo

# Analisis de simetria

#library(psych)
skew(horarios$FF)

[1] 1.591713

kurtosi(horarios$FF)

[1] 7.23501

La medida de asimetria y kurtosi terminan de validar lo que observamos en el histograma. La variable FF es asimetrica a derecha y tiene una mayor concentracion de valores muy cerca de la media de la distribución y muy lejos de la cola de la distribucion.

Necesitamos convertir las latitudes y longitudes a un valor decimal para poder graficarlo en el mapa con leaflet

#library(biogeo)

latitud <- c()
longitud <- c()

for(i in 1:nrow(estaciones)) {
     latitud[i] <- dms2dd(dd = estaciones[i, "LATITUD_GRADOS"], mm = estaciones[i, "LATITUD_MINUTOS"], ss = 0, ns = "S")
     longitud[i] <- dms2dd(dd = estaciones[i, "LONGITUD_GRADOS"], mm = estaciones[i, "LONGITUD_MINUTOS"], ss = 0, ns = "S")
}

estaciones['LATITUD'] <- latitud
estaciones['LONGITUD'] <- longitud

Unimos las dos tablas usando la variable NOMBRE como punto para combinar los datasets

data <- full_join(estaciones, horarios, by = c("NOMBRE" = "NOMBRE"))

glimpse(data)

Rows: 2,040
Columns: 18
$ NOMBRE           <chr> "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRANO II", "BASE BELGRA…
$ PROVINCIA        <chr> "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "ANTARTIDA", "…
$ LATITUD_GRADOS   <int> 77, 77, 77, 77, 77, 77, 77, 77, 62, 62, 62, 62, 62, 62, 62, 62, 63, 63, 63, 63, 63, 63, 63, …
$ LATITUD_MINUTOS  <int> 52, 52, 52, 52, 52, 52, 52, 52, 14, 14, 14, 14, 14, 14, 14, 14, 24, 24, 24, 24, 24, 24, 24, …
$ LONGITUD_GRADOS  <int> 34, 34, 34, 34, 34, 34, 34, 34, 58, 58, 58, 58, 58, 58, 58, 58, 56, 56, 56, 56, 56, 56, 56, …
$ LONGITUD_MINUTOS <int> 34, 34, 34, 34, 34, 34, 34, 34, 38, 38, 38, 38, 38, 38, 38, 38, 59, 59, 59, 59, 59, 59, 59, …
$ ALTURA           <int> 256, 256, 256, 256, 256, 256, 256, 256, 11, 11, 11, 11, 11, 11, 11, 11, 24, 24, 24, 24, 24, …
$ NRO              <int> 89034, 89034, 89034, 89034, 89034, 89034, 89034, 89034, 89053, 89053, 89053, 89053, 89053, 8…
$ NroOACI          <chr> "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYB", "SAYJ", "SAYJ", "SAYJ", "SAY…
$ LATITUD          <dbl> -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -77.86667, -62.…
$ LONGITUD         <dbl> -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -34.56667, -58.…
$ FECHA            <int> 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20042021, 20…
$ HORA             <int> 0, 3, 6, 9, 12, 15, 18, 21, 0, 3, 6, 9, 12, 15, 18, 21, 0, 3, 6, 9, 12, 15, 18, 21, 0, 1, 2,…
$ TEMP             <dbl> -20.0, -16.2, -16.9, -13.0, -15.1, -16.6, -15.3, -16.8, -0.1, -0.5, -0.9, -0.5, 0.3, 1.1, 1.…
$ HUM              <int> 70, 72, 80, 73, 86, 66, 77, 85, 88, 82, 77, 75, 84, 92, 99, 91, 43, 49, 67, 81, 67, 62, 77, …
$ PNM              <dbl> 975.7, 970.4, 966.7, 964.2, 961.4, 963.0, 964.8, 967.1, 992.9, 993.0, 992.9, 994.6, 995.3, 9…
$ DD               <int> 140, 140, 110, 140, 160, 160, 90, 110, 230, 230, 270, 290, 290, 320, 320, 320, 160, 200, 50,…
$ FF               <dbl> 13, 30, 50, 52, 74, 61, 44, 13, 15, 15, 26, 28, 41, 37, 50, 74, 20, 7, 13, 7, 11, 46, 56, 74…

Eliminamos las observaciones nulas en el dataset

data = na.omit(data)

Volvemos a graficar el grafico de correlacion incluyendo la variable altura.

corr <- cor(data[, c(7,13,14,15,16,17,18)], use = "complete.obs")
ggcorrplot(corr, type = "lower", outline.col = "black",
 lab=TRUE,
 ggtheme = ggplot2::theme_gray,
 colors = c("#6D9EC1", "white", "#E46726"))

Agrupamos los datos por nombre calculando el promedio y desvio del viento

data_agg = data %>%
  group_by(NOMBRE) %>%
  summarise(MEAN_VIENTO_KMH = mean(FF), SD_VIENTO_KMH = sd(FF), LONGITUD = unique(LONGITUD), LATITUD = unique(LATITUD), .groups = "keep")

Convertimos data_agg a data.frame ya que necesitamos este tipo de dato para poder trabajar

df_data_agg = data.frame(data_agg)

Transformamos df_data_agg en un archivo geográfico utilizando el código de proyección mercator

data_sf = st_as_sf(df_data_agg, coords = c("LONGITUD", "LATITUD"), crs = 4326)

as_Spatial(data_sf)

Error in as_Spatial(data_sf) : object 'data_sf' not found

Validamos la clase del nuevo dataframe

class(data_sf)

[1] "sf"         "data.frame"

##TODO: graficar el mapa de densidades con los puntos de las estaciones en data.

Grafico Interactivo

En el siguiente grafico de la republica argentina se observan en color azul las estaciones meteorologicas donde se realizaron las mediciones de la variables

#library(leaflet)
# Queremos que en el mapa se vea como etiqueta el nombre de la base meteorologica. Para eso aplicamos la siguiente funcion:
labs <- lapply(seq(nrow(data_sf)), function(i) {
  paste0( '<p>', data_sf[i, "NOMBRE"], '</p>' ) })

leaflet() %>%
  addTiles() %>%
  addCircles(data = data_sf, weight = 3, label = lapply(labs, htmltools::HTML))

NA

En el mapa observamos que hay puntos muy distantes de la argentina continental. Dado el proposito de este estudio, el cual es determinar la locacion geografica optima en base a la variable velocidad del viento, decidimos remover estas observaciones ya que no aporta informacion util y ademas agregan ruido a nuestro analisis.

Primero, vamos a borrar las estaciones que no estan en la argentina continental - Base Carlini - Base San Martin - Base Marambio - Base Esperanza - Base Orcadas

data_sf = data_sf[data_sf$NOMBRE != "BASE CARLINI (EX JUBANY)",]
data_sf = data_sf[data_sf$NOMBRE != "BASE SAN MARTIN",]
data_sf = data_sf[data_sf$NOMBRE != "BASE MARAMBIO",]
data_sf = data_sf[data_sf$NOMBRE != "BASE ESPERANZA",]
data_sf = data_sf[data_sf$NOMBRE != "BASE ORCADAS",]
data_sf = data_sf[data_sf$NOMBRE != "BASE BELGRANO II",]

Tambien vimos que El Bolson tiene las 9 observaciones en 0 para la variable FF y DD. Consideramos esto como un error en el instrumento de medicion, por lo cual vamos a eliminar a esa estacion del analisis.

data_sf = data_sf[data_sf$NOMBRE != "EL BOLSON AERO",]

Repetimos el plot para validar

labs <- lapply(seq(nrow(data_sf)), function(i) {
  paste0( '<p>', data_sf[i, "NOMBRE"], '</p>' ) })


leaflet() %>%
  addTiles() %>%
  addCircles(data = data_sf, weight = 3, label = lapply(labs, htmltools::HTML))

NA

Analisis Estadistico de los datos

Analisis univariado

describe(data_sf$MEAN_VIENTO_KMH)

hist(data_sf$MEAN_VIENTO_KMH)

boxplot(data_sf$MEAN_VIENTO_KMH)

hist(data_sf$SD_VIENTO_KMH)

boxplot(data_sf$SD_VIENTO_KMH)

#COMMENTS: sacar inliers (manual de buenas practicas compartido por silvia)

Analisis descriptivo de normalidad de los datos

hist(data_sf$MEAN_VIENTO_KMH)

boxplot(data_sf$MEAN_VIENTO_KMH)

qqnorm(data_sf$MEAN_VIENTO_KMH)
qqline(data_sf$MEAN_VIENTO_KMH, col=2)

#TODO: sacar outliers y puntos en la arg continental.

Hacemos el test de normalidad contra la variable MEAN_VIENTO_KHM

shapiro.test(data_sf$MEAN_VIENTO_KMH)

    Shapiro-Wilk normality test

data:  data_sf$MEAN_VIENTO_KMH
W = 0.9698, p-value = 0.03265

El p-value obtenido es menor a 0.05, lo cual indica que los datos que tenemos no son normales

Eliminamos outliers

Inliers

knea <- knearneigh(data_sf)
neib <- knn2nb(knea)


listw <- nb2listw(neib)

globalMoran <- moran.test(data_sf$MEAN_VIENTO_KMH, listw)
globalMoran

    Moran I test under randomisation

data:  data_sf$MEAN_VIENTO_KMH  
weights: listw    

Moran I statistic standard deviate = 4.0432, p-value = 2.636e-05
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.52033029       -0.01111111        0.01727672 

moran <- moran.plot(data_sf$MEAN_VIENTO_KMH, listw = listw)

data_sf[10,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -65.76667 ymin: -28.6 xmax: -65.76667 ymax: -28.6
Geodetic CRS:  WGS 84
           NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                geometry
16 CATAMARCA AERO        46.28571      11.51397 POINT (-65.76667 -28.6)

data_sf[56,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -62.38333 ymin: -37.6 xmax: -62.38333 ymax: -37.6
Geodetic CRS:  WGS 84
       NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                geometry
63 PIGUE AERO        38.23077      6.905962 POINT (-62.38333 -37.6)

data_sf[87,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -61.96667 ymin: -33.66667 xmax: -61.96667 ymax: -33.66667
Geodetic CRS:  WGS 84
               NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                    geometry
94 VENADO TUERTO AERO              32      11.09955 POINT (-61.96667 -33.66667)

data_sf[21,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -61.88333 ymin: -37.43333 xmax: -61.88333 ymax: -37.43333
Geodetic CRS:  WGS 84
                NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                    geometry
27 CORONEL SUAREZ AERO        27.33333      5.278889 POINT (-61.88333 -37.43333)

data_sf[73,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -66.35 ymin: -33.26667 xmax: -66.35 ymax: -33.26667
Geodetic CRS:  WGS 84
          NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                 geometry
80 SAN LUIS AERO              17      7.901568 POINT (-66.35 -33.26667)

data_sf[51,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -64.31667 ymin: -23.15 xmax: -64.31667 ymax: -23.15
Geodetic CRS:  WGS 84
      NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                 geometry
58 ORAN AERO        1.222222      2.981424 POINT (-64.31667 -23.15)

data_sf[55,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -71.01667 ymin: -46.51667 xmax: -71.01667 ymax: -46.51667
Geodetic CRS:  WGS 84
               NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                    geometry
62 PERITO MORENO AERO             1.1      3.478505 POINT (-71.01667 -46.51667)

data_sf = data_sf[data_sf$NOMBRE != "CERES AERO",]
data_sf = data_sf[data_sf$NOMBRE != "SAN LUIS AERO",]
data_sf = data_sf[data_sf$NOMBRE != "PERITO MORENO AERO",]
data_sf = data_sf[data_sf$NOMBRE != "ORAN AERO",]
data_sf = data_sf[data_sf$NOMBRE != "CORONEL SUAREZ AERO",]
data_sf = data_sf[data_sf$NOMBRE != "PIGUE AERO",]
data_sf = data_sf[data_sf$NOMBRE != "CATAMARCA AERO",]

Inliers: verificamos los resultados dsp de eliminar inliers

# Creamos una lista de vecinos
knea <- knearneigh(data_sf)
neib <- knn2nb(knea)
listw <- nb2listw(neib)

# Hacemos el test de moran 
globalMoran <- moran.test(data_sf$MEAN_VIENTO_KMH, listw)
globalMoran

    Moran I test under randomisation

data:  data_sf$MEAN_VIENTO_KMH  
weights: listw    

Moran I statistic standard deviate = 4.8984, p-value = 4.83e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.66907198       -0.01204819        0.01933449 

moran <- moran.plot(data_sf$MEAN_VIENTO_KMH, listw = listw)

data_sf[81,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -63.01667 ymin: -40.85 xmax: -63.01667 ymax: -40.85
Geodetic CRS:  WGS 84
        NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                 geometry
95 VIEDMA AERO            16.5      5.890302 POINT (-63.01667 -40.85)

data_sf[60,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -69.28333 ymin: -51.61667 xmax: -69.28333 ymax: -51.61667
Geodetic CRS:  WGS 84
              NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                    geometry
73 RIO GALLEGOS AERO        25.91667      6.907944 POINT (-69.28333 -51.61667)

data_sf[80,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -61.96667 ymin: -33.66667 xmax: -61.96667 ymax: -33.66667
Geodetic CRS:  WGS 84
               NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                    geometry
94 VENADO TUERTO AERO              32      11.09955 POINT (-61.96667 -33.66667)

data_sf[19,]

Simple feature collection with 1 feature and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -58.76667 ymin: -27.45 xmax: -58.76667 ymax: -27.45
Geodetic CRS:  WGS 84
            NOMBRE MEAN_VIENTO_KMH SD_VIENTO_KMH                 geometry
28 CORRIENTES AERO        9.916667      5.115336 POINT (-58.76667 -27.45)

shapiro.test(data_sf$MEAN_VIENTO_KMH)

    Shapiro-Wilk normality test

data:  data_sf$MEAN_VIENTO_KMH
W = 0.96946, p-value = 0.04325

data_sf = data_sf[data_sf$NOMBRE != "VIEDMA AERO",]
data_sf = data_sf[data_sf$NOMBRE != "RIO GALLEGOS AERO",]
data_sf = data_sf[data_sf$NOMBRE != "VENADO TUERTO AERO",]
data_sf = data_sf[data_sf$NOMBRE != "CORRIENTES AERO",]

# Creamos una lista de vecinos
knea <- knearneigh(data_sf)
neib <- knn2nb(knea)
listw <- nb2listw(neib)

# Hacemos el test de moran 
globalMoran <- moran.test(data_sf$MEAN_VIENTO_KMH, listw)
globalMoran

    Moran I test under randomisation

data:  data_sf$MEAN_VIENTO_KMH  
weights: listw    

Moran I statistic standard deviate = 4.8722, p-value = 5.518e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.67737208       -0.01265823        0.02005783 

moran <- moran.plot(data_sf$MEAN_VIENTO_KMH, listw = listw)

plot(data_sf$MEAN_VIENTO_KMH)
plot(log(data_sf$MEAN_VIENTO_KMH))

qqnorm(log(data_sf$MEAN_VIENTO_KMH))
qqline(log(data_sf$MEAN_VIENTO_KMH), col=2)

shapiro.test(data_sf$MEAN_VIENTO_KMH)

    Shapiro-Wilk normality test

data:  data_sf$MEAN_VIENTO_KMH
W = 0.97049, p-value = 0.06096

Variograma

#coordinates(data_sf) = ~LATITUD+LONGITUD
#class(data_sf)

plot(df_data_agg)
plot(data_agg[, c(4,5)])

v <- variogram(MEAN_VIENTO_KMH~1, data_sf)
plot(v)

v <- variogram(MEAN_VIENTO_KMH~1, data_sf)
plot(v)

vt_exp = fit.variogram(v, vgm(125, "Exp", 30, 5))
vt_exp

plot(v , vt_exp)

vt_mat = fit.variogram(v, vgm(125, "Mat", 30, 5))
plot(v , vt_mat)

vt_exc = fit.variogram(v, vgm(125, "Exc", 30, 5))

No convergence after 200 iterations: try different initial values?

plot(v , vt_exc)

Kriging

attr(vt_exp, 'SSErr')

[1] 0.1189922

attr(vt_mat, 'SSErr')

[1] 0.1189922

attr(vt_exc, 'SSErr')

[1] 11.92923

attr(vt_bes, 'SSErr')

[1] 0.1344978

ko1 <- krige(MEAN_VIENTO_KMH~1, data_sf, grilla, model = vt_exp, nmax=20)

spplot(ko1["var1.pred"])
spplot(ko1["var1.var"])